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In his famous 1981 talk, Feynman proposed that unlike classical computers, which would presumably expe- 
rience an exponential slowdown when simulating quantum phenomena, a universal quantum simulator would 
not. An ideal quantum simulator would be error resistant, easily controllable, and built using existing technol- 
ogy. Moving away from gate-model and projective measurement based implementations of quantum computing 
may offer a less resource-intensive, and consequently a more feasible solution. Here we consider an adiabatic 
quantum simulator which simulates the ground state properties of sparse Hamiltonians consisting of one- and 
two-body interaction terms, using sparse Hamiltonians with at most three-body interactions. Properties of such 
Hamiltonians can be well approximated with Hamiltonians containing only two-local terms. The register hold- 
ing the simulated ground state is brought adiabatically into interaction with a probe qubit, followed by a single 
diabatic gate operation on the probe which then undergoes free evolution until measured. This allows one to 
recover e.g. the ground state energy of the Hamiltonian being simulated. Given a ground state, this scheme can 
be used to verify the QMA-complete problem LOCAL HAMILTONIAN, and is therefore likely more powerful 
than classical computing. 

I. INTRODUCTION 

Computer simulation of quantum mechanical systems is an indispensable tool in all physical sciences dealing with 
nanoscale phenomena. Except for specific and rare casesi, classical computers have not been able to efficiently 
simulate quantum systems, as in all known techniques at least one of the computational resources required to perform 
the simulation scales exponentially with the size of the system being simulated. 

Numerous classical approximative methods, such as density functional theory (DFT 2 ) and quantum Monte Carlo 
(QMC 3 ) have been developed to address various aspects of the efficiency problem, but no known polynomially scaling 
methods are universally applicable. Each suffers from particular deficiencies such as the fermionic sign problem of 
QMC or the approximate exchange-correlation functionals of DFT. Quantum computers on the other hand, as conjec- 
tured by Feynman 4 , may be used to simulate other quantum mechanical systems efficiently. Feynman's conjecture was 
subsequently proven and expanded leading to the rapidly growing area of study known as quantum simulation^—. 

Quantum simulation is expected to be able to produce classically unattainable results in feasible run times, using 
only a modest number of quantum bits. For example, calculating the ground state energy of the water molecule to 
the level of precision necessary for experimental predictions (rs 1 kcal/mol) — a problem barely solvable on current 
supercomputers— — would require roughly 128 coherent quantum bits on a gate-model quantum computer 9 . However, 
this calculation would also require on the order of billions of quantum gates and thus quantum error correction would 
be essential' 4 . The large number of gates and qubits needed for practical applications of quantum simulation has lead 
current experimental work to fine tune small system s 11 - 15 ' 16 . 

Theoretical quantum simulation falls into two main categories: dynamic evolution and static properties. Both cat- 
egories rely heavily on the Trotter decomposition to handle non-commuting terms in the Hamiltonian when mim- 
icking the unitary time propagator of the system to be simulated. To approximate evolution under a Hamilto- 
nian H = Ylt=i Hi consisting of k non-commuting but local terms {Hi\\ =1 , one applies the sequence of unitary 
gates { e -*- H 'i t / n }|_ 1 a total of n times. As number of repetitions n tends to infinity the error caused by the non- 
commutativity vanishes and the approximation converges to the exact result 5 . If each time slice requires a constant 
number of gates independent of the parameter t/n, then reducing the error by repeating the sequence n times can 
become prohibitively expensive for high accuracy applications. 

Constructing a practical method of quantum simulation is a significant challenge (see in particular the recent work 14 
which analyzes fault tolerance in a gate-model quantum simulation algorithm). Moving away from the gate model 
may offer a less resource-intensive, and consequently a more feasible solution. This letter addresses the problem by 
presenting a hybrid model of quantum simulation, consisting of an adiabatically controlled simulation register coupled 
to a single gate-model readout qubit. Our scheme can simulate the ground state properties of arbitrary spin graph 
Hamiltonians. Furthermore, it allows the measurement of the expectation value of any constant fc-local observable 
using a k + 1-local measurement Hamiltonian. Even in the presence of considerable noise at least some information 
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can be extracted. 

In most quantum computing architectures, the natural interactions are 2-body. However, under certain conditions k- 
body interactions can be well approximated using techniques such as perturbative Hamiltonian gadgets 17-19 or average 
Hamiltonian methods 20 . In particular, reference 21 considered the mapping between a given n-qubit target Hamiltonian 
with A: -body interactions onto a simulator Hamiltonian with two-body interactions. 

Structure of this paper: We will continue by giving an overview of the simulator design, including algorithm 
initialization, adiabatic evolution and measurement. One of the key motivations for the proposed model is its resistance 
to certain types of noise which is outlined in Section [HI] and we conclude in Section [IV] The supporting material in 
the Appendix further details the specifics of our method (including the full numeric simulation of the technique in the 
presence of noise). 

H. SIMULATOR OVERVIEW 

We consider the simulation of systems represented by finite collections of spins acted on by a time-independent 
Hamiltonian described by a graph G — (V, E) — e.g. the Heisenberg and Ising models. Each graph vertex v € V 
corresponds to a spin acted on by a local Hamiltonian L v , and each edge e g E to a pairwise interaction K e between 
the involved vertices. The Hamiltonian Ht we wish to simulate is given by 

The simulator consists of an adiabatically controlled simulation register S with Hamiltonian Hs, and a probe regis- 
ter P which will be acted on by gate operations and measured projectively. We will engineer the probe P such that it 
behaves as a controllable two-level system with orthonormal basis {\po), \pi)}- Without loss of generality, the probe 
Hamiltonian can be expressed as Hp — 8\pi)(pi\, where 6 is the spectral gap between the probe's ground and first 
excited states. 

Initialization: We will first set the simulation register Hamiltonian Hs to Hgj and prepare S and P in then- 
respective ground states. The Hamiltonian Hs,i has a simple ground state which can be (i) computed classically 
and (ii) prepared experimentally in polynomial time — such as a classical approximation to the ground state of the 
simulated system. The simulator Hamiltonian is thus initially given by 

#o = H S j ® t P + t s ® Hp. (2) 

Adiabatic evolution: According to the adiabatic theorem 22 , a quantum system prepared in an energy eigenstate will 
remain near the corresponding instantaneous eigenstate of the time-dependent Hamiltonian governing the evolution 
if there are no level crossings and the Hamiltonian varies slowly enough. By adjusting the simulator parameters, we 
adiabatically change Hs from Hsj to Hs,t, the fully interacting Hamiltonian of the system to be simulated. 

Let us denote the ground state of Hs.t as \sq). At the end of a successful adiabatic evolution P is still in its ground 
state \po), and S is in (a good approximation to) the ground state of the simulated system, \so). Hence the simulator is 
now in the ground state |g) = |so) ® \Po) of its instantaneous Hamiltonian 

Hi = H S .T ® lp + Is <8> H P . (3) 

The computational complexity of preparing ground states of quantum systems has been studied^—. It is possible 
to prepare a desired ground state efficiently provided that the gap between the ground and excited states is sufficiently 
large 22 . This depends on the adiabatic path taken. In general, finding the ground state energy of a Hamiltonian, 
even when restricted to certain simple models, is known to be complete for QMA, the quantum analogue of the class 
jj enC g 5 it i s considered unlikely that even a quantum computer would be able to do this in polynomial time. 
In fact there are physical systems such as spin glasses in nature which may never settle in their ground state. However, 
it is hoped (but not proven) that a host of realistic systems will satisfy the property of having a large energy gap on 
even simple linear adiabatic paths and thus be amenable to quantum simulation algorithms which rely on adiabatic 
state preparation. 

Measurement: The measurement procedure begins by bringing S and P adiabatically into interaction. The simulator 
Hamiltonian becomes 



H 2 = Hi + A®\ Pl )(p 1 \ , 

H S p 



(4) 
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where the operator A corresponds to an observable of the simulated system that commutes with Hs.t- Assuming 
that A can be decomposed into a sum of two-local operators, the interaction term H$p involves three-local interac- 
tions. These terms can be implemented using either Hamiltonian gadget techniques or average Hamiltonians (see the 
supporting material for details). 

Let us use ao := (sq\A\so) to denote the expectation value we wish to measure. If ao + 5 > 0, then |g) is also 
the ground state of H2 (see the supporting material for proof). At time t = we apply a Hadamard gate to the 
measurement probe which puts it into a superposition of its two lowest states. This is no longer an eigenstate of H2, 
and the system will evolve as 

|^(t)) = ^K>®(|po)+e- iwt |p 1 )), (5) 

where w := (a + S)/H. We have thus encoded the quantity a , a property of the ground state of Hst> into the time 
dependence of the probe P. 

After a time t, we again apply a Hadamard pulse to the probe, resulting in the state 

|V»(*)> = l*o) ® (cos (wt/2) bo) + i sin (wt/2) \pi}) , (6) 
and then measure the probe. The probability of finding it in the state \po) is 

^b(i) = 5 0- + cos(urf)) = cos 2 (wt/2) . (7) 

If we have non-demolition measurements (see e.g.—) at our disposal, then measuring the probe does not disturb the 
state of the simulator which can be reused for another measurement. 

One repeats the measurement with different values of t until sufficient statistics have been accumulated to recon- 
struct uj and hence ao — this is reminiscent of Ramsey spectroscopy 24 and hence should seem natural to experimen- 
talists. In essence, we have performed Kitaev's phase estimation algorithm 25 , using the interaction Hamiltonian Hsp 
instead of the controlled unitary. 

If the ground state subspace of Hs,t is degenerate and overlaps more than one eigenspace of A, or the simulation 
register S has excitations to higher energy states at the beginning of the measurement phase, the probability of finding 
the probe in the state |po) is given by a superposition of harmonic modes. For example, for a thermalized state with 
inverse temperature f3, we obtain 




-L^^e-^cosK,*)) 
x v ki / 



(8) 



where the first summation index runs over the energy eigenstates and the second over the eigenstates of A in which 
energy has value Ek, and Wk,i — (a,k,t + S)/H. 



III. EFFECTS OF NOISE 



To assess the effects of noise on the simulation method, we performed a numerical simulation of the simplest 
nontrivial implementation of the hybrid simulator, consisting of two simulator qubits and one probe qubit, with a 
randomly chosen Ht- Each qubit was coupled to its own bosonic heat bath with an Ohmic spectral density using the 
Born-Markov approximation 26 . The qubit-bath couplings were chosen such that the resulting single-qubit decoherence 
times T\ and T2 are compatible with recent superconducting flux qubit experiments with fully tunable couplings (see 
for example 27 ). The noise model is described further in the supporting material. The observable A was chosen to be 
Hs,t, the simulated Hamiltonian itself, which means that the ground state energy was being measured. 

We simulated measurements on 40 evenly distributed values of the time delay t. At each t.- L we performed 50 
measurements, and averaged the results to get an estimate of Po(tj). An exponentially decaying scaled cosine function 
was then fitted to this set of data points to obtain an estimate for w and thus for sq. The fit was done using the MATLAB 
LSQNONLIN algorithm, guided only by fixed order-of-magnitude initial guesses for the parameter values. 

The results of the simulation are presented in Fig. Q] The noise, together with the slight nonadiabaticity of the 
evolution, cause excitations out of the ground state which result in a signal consisting of multiple harmonic modes. 
However, the ground state mode still dominates and with a realistic level of noise and relatively modest statistics 
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Ramsey measurement of the probe 




FIG. 1: Measurement procedure under Markovian noise. The continuous curve represents the probability of finding the probe in 
the state |po), the circles averaged measurement results and the dotted line a least squares fit to them, fkuo — h ■ 25 MHz is the 
energy scale of the Hamiltonians involved. 

we are able to reconstruct u to a relative precision of better than 0.01. This includes both the uncertainty due to 
finite statistics, and the errors introduced by the environment through decoherence and the Lamb-Stark shift. In an 
experimental implementation there may also be other noise sources not considered here related to the measurement 
process itself, as well as systematic (hardware) errors in the quantum processor (qubits, couplers etc.). 

IV. CONCLUSION 

The presented simulation scheme significantly differs from gate-model methods in that it does not require millions 
of coherent gate operations. Instead, an adiabatic control sequence is used, which requires much less complicated 
control hardware and is likely to be more robust against errors. A simple experimental implementation could be 
feasible with present-day technology. A direct comparison to a gate-model simulator using phase estimation to extract 
the result is difficult because such a design would have to use quantum error correction to be useful at the noise levels 
we have considered. At the measurement stage we do require single-qubit gate operations and measurements, but 
only on the probe qubit. Furthermore, these operations are relatively simple to implement compared to a full Trotter 
decomposition of the simulated Hamiltonian. It is well known that fault tolerance in the adiabatic model depends on 
the energy difference between the Hamiltonian's ground and first excited states. Several methods for the adiabatic 
model have been proposed which could in principle be used to increase the fault tolerance of our protocol^—. 

In order to simulate a system of n qubits with a 2-local Hamiltonian described by the graph G, ideally our scheme 
requires one probe qubit for the readout and n simulation qubits. Additionally, if Hamiltonian gadgets are used 
to implement 3-local interactions, one ancilla qubit is required for each two-local term in the simulated observable 
(represented by an edge in the corresponding graph). In practice the number of ancillas may be slightly higher if more 
involved types of gadgets are used to implement the 3-local interactions. The total number of qubits required thus 
scales as 0(n) for sparse graphs and 0(n 2 ) for maximally connected ones. 
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Appendix A: Supporting material 

We have created a MATLAB-based numerical simulation of a simple case of the presented hybrid simulation method 
under Markovian noise. The full source code of the simulation will be made available online through EPAPS. 



1. The Gadget 

As a proof of concept, our simulator uses 2-body interactions to simulate a Hamiltonian with a 3-body interaction. 
This Gadget Hamiltonian construction was proposed in 17 and has since been used and extended by other a 18,19 . These 
papers contain further background details including notation. 

We label the working qubits 1-3. The effective 3-body interaction is well approximated on qubits 1, 2 and 3 
(within e) in a low energy subspace — constructed by adding a penalty Hamiltonian to an ancillary (mediator) qubit to. 
The mediator qubit to doubles the size of the state space and the penalty Hamiltonian splits the Hilbert space into low 
and high energy subspaces, separated by an energy gap A (which is inversely proportional to a polynomial in e — in 
our case reduced to A > e~ 3 ). The details of the gadget we develop for use in this work follow. 
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We will apply a Hamiltonian H p to the mediator qubit m as well as a Hamiltonian V to qubits 1-3 and m. The 
Hamiltonian H p + V has a ground state that is e-close to the desired operator J A\ ® A2 <S> A3 ® |0) (0|. 

H p :=A|l)(l| m (Al) 



V(J, A(e)) := y + A 1/,3 |0)(0| m - A 1 / 3 ^ ® A 2 + - A x ) ®a x + JA 3 ® (1 - A 2 / 3 |l)(l| m ) (A2) 

where y is some Hamiltonian already acting on qubits 1 — 3 as well as a possible larger Hilbert space. Note the gadget 
above assumes Af = 1, Vi = 1, 2, 3. 

The so called self-energy expansion (IA3b under appropriate conditions is known to provide a series approximation 
to the low-lying eigenspace of an operator. To verify that the Hamiltonian H p + V gives the desired approximation, 
one relies on expansion of the self energy to A th order (here the higher order terms give rise to effective interactions 
greater than second order): 

m.) - wm + warn + n|ma +P (gp) , (A „ 

where the operator is written in the {|0), |1)} basis of m. One considers the range \z\ < 2\J\ + e and notes that 
||S_(z) - J Ax ® A 2 ® A 3 1| = 0(e) and applies the Gadget Theorem^. 

Before concluding this section, we note that care must be taken when adiabatically evolving gadgets. Reference^ 
contains a proof that the linear path Hamiltonian is universal for adiabatic quantum computation. Universality (and 
hence a non-exponentially contracting gap) remains when the locality of the construction is reduced using perturbative 
gadgets^. 



2. Average Hamiltonian Method 



An alternate method of generating the special Hamiltonians we require is to make use of the time average of a 
series of simpler generating Hamiltonians 20 . It has long been known that by regularly switching between a set of fixed 
Hamiltonians {Hi}, it is possible to approximate time evolution under any other Hamiltonian H, provided that H is 
contained within the algebra generated by {Hi}. This fact lies at the heart of both average Hamiltonian theory and 
geometric control theory. Over the years many methods have been developed for making the approximations accurate 
to high order, however here we will focus only on approximations correct to first order. 

In order to construct an average Hamiltonian we will make use of a first order approximation to the Baker-Campbell- 
Hausdorff formula: 

log(e A e s ) ps A + B + ^[A,B]. (A4) 

From this we obtain formulae for approximating the exponential of both the sum and the Lie bracket of A and B to 
first order. 

e A+B « e A e B , (A5) 
e [A,B] w e A e B e ~A e -B_ (A6) 



These equations can be related to time evolution under some Hamiltonians, Hi, by replacing A and B with operators 
of the form iHit/h. Clearly by applying these rules recursively, it is possible to generate any Hamiltonian in the Lie 
algebra generated by the initial set of static Hamiltonians. We note that the combination of any pairwise entangling 
Hamiltonian, such as an Ising Hamiltonian together with tunable local Z and X fields is sufficient to generate the full 
Lie algebra su(2 N ) for an N qubit system, and so can be used to approximate an arbitrary Hamiltonian. 

Although the time-varying Hamiltonian means that the system does not have a ground state in the normal sense, 
the average energy of the system is minimized when the system is within 0(t) of the ground state of the target 
average Hamiltonian. As a result of this, if the time scale for switching between Hamiltonians is small compared to 
the time scale for the adiabatic evolution of the system, the system will behave as if it was experiencing the average 
Hamiltonian. 
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TABLE I: Noise model parameters 

Ti ~ N(1.0, 0.1) fis 

T<f, ~ iV(1.3, 0.1) ijs 

T2 1 =\T- 1 +T- 1 

T 20 mK 

LO c 20 Wo 

cjo 2tt • 25 MHz 



3. Noise model 



The noise model used in our MATLAB simulation consists of a separate bosonic heat bath coupled to each of the 
qubits. The baths have Ohmic spectral densities, 

J(oj) = h 2 uje{uj)e{uj c - uj), (A7) 

where the cutoff uj c was chosen to be above every transition frequency in the system, and the baths are assumed to be 
uncorrelated. Each bath is coupled to its qubit through an interaction operator of the form 

D = X(cos(a)a z + sin(a)cr x ). (A8) 

Using the Born-Markov approximation we obtain an evolution equation which is of the Lindblad form 26 . Denoting 
the level splitting of a qubit by HA, we obtain the following uncoupled single-qubit decoherence times: 

Tf 1 = A 2 sm 2 (a)27rAcotli(^A/2), (A9) 
T 2 _1 = iT-f 1 + A 2 cos 2 (a)47r/(^), (AlO) 

where /3 = ■^ T . Given T U T 2 ,T and A, we can solve the bath coupling parameters A and a separately for each qubit, 
and then use the same noise model in the fully interacting case. The values used for the parameters are presented 
in Table Q] The single-qubit decoherence times T± and T2 were chosen to be compatible with recent coupled-qubit 
experiments with fully tunable couplings such as 27 . The bath temperature T and the energy scale huio for of the 
effective Hamiltonian to be simulated were likewise chosen to match those found in typical superconducting qubit 
experiments. To simulate manufacturing uncertainties, the actual values of the T\ and parameters for the individual 
qubits were drawn from a Gaussian distribution with a small standard deviation. 

4. Measurement 

The pre-measurement system Hamiltonian is 

Hi = H s ,t ® 1 + 1 <8> Hp. (All) 

The operators Hs and Hp can be expanded in terms of their eigenvalues and eigenstates using the (possibly degener- 
ate) spectral decomposition 

H s = ^2 s k \s k j)(s k ,j\, s <si<... (A12) 

kj 

and correspondingly for Hp. 

Let the states of systems S and P begin in their ground state subspaces, the full normalized state |g) of the simulator 
belonging to the ground state subspace of the non-interacting Hamiltonian H\\ 

|g) G span{|s 0!fc ) ® \p ,i}} k ,i- (A13) 

S and P are brought adiabatically into interaction with each other. The Hamiltonian becomes 

H 2 =Hi+A®(t P -Tl Po ), (A14) 
s , ' 

Hsp 
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where the operator A corresponds to an observable of the simulated system that commutes with H$, and 

n po = ^2 \Po,m)(P0,m\ (A15) 
m 

is the projector to the ground state subspace of Hp. Because A and H$ commute, they have shared eigenstates: 

Hs\sk,j) = s k \sk,j), (A16) 
A\s htj ) = ak,j\sk,j)- (A17) 

a. A ground state lemma 

We will now show that Hamiltonians Hi and H2 have the same ground state subspace given that a m i n + pi — pa > 
where a m ; n is A's lowest eigenvalue. 

Lemma 1. Let Hi and H2 be the finite dimensional Hamiltonians defined previously in ( IA1 11 1 and ( IA141 >. and a m i„ + 
Pi — Po > 0. Now, iff Hi\*) — \\k), where A is the smallest eigenvalue of Hi, then H2I*) = l|*)> where k is the 
smallest eigenvalue of H2- 

Proof. Firstly, we have 

(g\Hi\g) = s +po, 
(g\H SP \g) = 0. 

Now, expanding an arbitrary normalized state \(j>) in the eigenbases of H$ and Hp, 

\4>) = y ] c xyij\s x ,i)\Py.j) , 



we get 



{(f>\Hl\4>) = ^ \Cx yl j\ 2 (s x +Py) 

xyij 

> ^ \ c xyij\ 2 ( s o +Po) = s o +Pa, 

xyij 

(4>\H SP \4>) = (0| Y c xyij A\s Xti )(lp - Ilp )\Py,j) 

xyij 

= (^1 y y c xyijQ-x,i\S X ,i) \Py,j) 
xij y>l 

xij y>l 

> a m in X! \ Cx vi]\ 2 - if a min > 0. 

xij y>l 

Hence, if A's lowest eigenvalue a m i n > then \g) is the ground state of H2 = Hi + Hsp as well. If A is not 
nonnegative, we can perform the transformation 

H's = Hs + Clminl, 

Hp = Hp — a m i n n po , 

A = A — amin/. 

This leaves H2 invariant, but makes A' nonnegative. As long as a m i„ + pi — po > then {\po y m )} m still span the 
ground state subspace of H' p and the above analysis remains valid. □ 
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b. Ground state degeneracy and excitations 

If the ground state subspace of Hs,t is degenerate and overlaps more than one eigenspace of A, or the simulation 
register S has excitations to higher energy states at the beginning of the measurement phase, we need a more involved 
analysis of the measurement procedure. Assume the pre-measurement-phase state of the simulator is given by 

Pa = ( £ c ki mn \skd){s m ,n\ j ® y 1 ° J • (A18) 

\klmn / \ / 

All three terms in H 2 commute given that [H s , A] = 0. Hence 

e qH 2 = ( e q Hs i p )(i s ® e qHp ){e qHsp ) 

= (e qHs ® e qHp )(e qA ® (t P - Il po ) + t s ®U po ). (A19) 
As a result of the measurement procedure right before the actual measurement the state is given by 

pi(t) = Upe- ltH2 / R Upp Upe itH2 / h Up 



\klmn 



= £ (c„ mB e-^-™ >' V«N-™,»l) ® ^Hp ( e _ itWM [ 2fl _ ^I^i^/i) ^, (A20) 

k I m n \ / 

where C0k,i = (a,k,l + S)/h and Hp is the Hadamard gate operating on the probe. Projecting the probe to the ground 
state subspace, we get H po pi(t)U pa = 

= (cki m ne- lt{Sk - s ™ )/h \ Sk j)( Sm , n \) 1(1 + e -*(°M-«m.»)/" + 2cos( Wm ,„t)(2a - 1))|0)(0| (A21) 
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Thus the probability of finding the probe in the ground state subspace is given by a superposition of harmonic modes: 

P (t) = TT(U PoPl (t)U Po ) = 1(1 + £ c fcIfc , cos(w fc ,,t)(2a - 1)). (A22) 



